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ABSTRACT 


For the first time a two dimensional flow model was' 
developed for a shell and tube heat exchanger without phase 
change. For the shell side fluid, the conservation equations 
of mass and momen-tum were modified using the concepts of 
volume porosity and distributed resistance to account for 
the presence of a tube bundle. The conservation equations 
were solved using the Finite Element Method to obtain the 
velocity field inside a shell and tube heat exchanger. The 
computational method used is capable of handling a finer 
mesh for a large calculation domain and it gives a faster 
convergence. It also makes it possible to save a lot of 
computer core memory and computational time. The results 
obtained were encouraging and showed an expected behavior. 

The flexibility of the computer code developed 
makes it useful to model flows in many other types of heat 
exchangers , 



CHAPTER-I 


INTRODUCTION 


1 .1 MOTIVATION : 

Fluid -dynamics theory, for reasons of history and 
fashion, has paid little attention to the phenomena which 
occiAT in industrial equipment such as heat exchangers, 
cooling towers etc. Even now, when the digital computers 
have removed the obstacles of numerical computation, the 
prized achievements of applied mathematicians concern more 
often the supersonic flow of air about a missile, or the 
development of an eddy behind a cylinder, than any phenomenon 
of interest to a process engineer. 

As a consequence, the designer of process equipment 
must usually base predictions of its performance on assump- 
tions about the flow patterns, departures from these 
assumptioris are then allowed for by empirically -derived 
correction formulae . Thus a single pass shell and tube heat 
exchanger may be designed by reference to the formxila for 
the performance of an ideal counter flow heat exchanger, 
modified by some performance degradation factor which takes 
account of departures from ideal flow pattern. 



2 


Unfortunately, although tliis practice can serve 
adequately for interpolation between equipment items of 
common type with only slight variations of geometry, it fails 
completely where the benefit of new geometrical configurations 
is to be explored; optimum configurations must, therefore, 
be determined ejqjerimentally as a rule . None tlie less doing 
an experiment using a model is not an easy task. Besides, 
it is a time consuming and costly affair compared with its 
equivalent computational work. 


The present work is an attempt to introduce a new 
fashion in flow modeling of process equipment, by applying 
to a shell and tube heat exchanger a numerical procedure 
for calculating the fluid -flow field. 

The shell side flow is considered as a two dimensional 
flow of liquid through a porous medium which consists of 
tube bundles . The flow-field on shell side is calculated 
using finite -element method. The tube side flow is assumed 
to be one -dimensional , Both the liquids on shell and tube 
side remain in liquid state throughout the space inside the 
heat exchanger. The vdiole analysis is done for steady -state 
operation of the heat exchanger. 
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1 .2 ORGANISATION OF THESIS : 

The organisation of thesis is a simple one . The 
thesis is divided into four chapters excluding the 
present one . 

Chapter two is a literature smvey related to the 
present problem and its relation to the present work 
undertaken. The objectives of the present work are also 
listed in this chapter. 

Chapter three describes the modeling techniques 
used for finding out the velocity fields on shell 
side of the heat exchanger. The governing differential 
equations and the formulation are discussed here. 

Chapter four describes the solution procedure and 
the programme for the calculation of velocity field inside 
the heat exchanger. The computational method employed for 
the solution of algebraic equations is also given in this 
chapter . 

Chapter five gives the discussion of results and 
the conclusions . The limitations of the present work and 
the scope for future work are also discussed in this 
chapter. This chapter is followed by a list of references. 
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CHAPTER- II 
LITERATURE SURVEY 

2.1 REVIEW OF THE EXISTING MODELS : 

A few simplified models are available for the 
steady -state analysis of shell and tube heat exchangers. 
Almost all of them are developed for the analysis of 
steam generators. 

The fundamental conservation principles are applied 
in all the models developed so far , Ihe differences in the 
various formulations available lie in the way the conserva- 
tion equations are written, the assumptions and the approxi- 
mations in the equations , The models available are either 
one -dimensional or three-dimensional and use Finite - 
Difference technique in most of the cases, for the Solution 
of the problem. 

The idea of using distributed resistances to 
simulate the presence of heat transfer tubes and baffle 
plates on the shell side of a heat exchanger was first 
introduced by Patankar and Spalding /~’1974__7, They had 
assumed that the space inside a heat exchanger was uniformly 
filled with fluid, throughout which, however a resistance 
to fluid motion was distributed on a fine scale. They did 
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not calculate tute-side flow distribution. The model was 
developed using continuum approach by finite -difference 
scheme. Diffusion and shear-forces were not taken into 
account while simulating the model. Also the authors 
assumed constant density, throughout the calculation domain. 

More recently AbuRomia et al, applied the 

distributed resistance concept to obtain the flow field 
between the typical tube support spans of the Clinch River 
Breeder Reactor plant (CRBRP) Intermediate Heat Exchangers. 

Patankar et. al, developed a computer code 

which was a "three -dimensional finite -difference analysis 
for predicting the local thermal /hydraulics of nuclear 
ojnce -through generators. The tube bundle and tube -support 
plates were modeled as a porous media with distributed 
resistance to flow. They used quasi -continuimi approach and 
solved a steady state, fully elliptic set of governing 
eq-uations for the simulation, 

A complicated analysis of Shell and Tube Heat 
Exchangers using finite -difference method was given by 
Sha et. al._^l982_7. They also used continuum approach to 
solve Navier -Stoke ' s equations which were modified to 
represent the phenomenon in a heat exchanger, in cylindrical 
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co-ordinates. They considered volume porosity as well as 
surface permeability to model the flow. Since, in general, 
the tube bundle porosity and permeability are anisotropic, 
they used a combination of porosity, permeability and 
distributed resistance for the analysis of flow-field. 
Sirrface permeability was defined as the fraction of open 
projected flow area in the direction of flow component in 
the control volume. They had also developed a one-equation 
turbulence model for a tube bundle based on transport 
equation for the turbulent kinetic energy. The model took 
into account the effect of diffusion and shear forces. 

One dimensional analysis for steam generator 
transients was given by Bhatnagar /”1983__7 using finite 
difference technique, 

2.2 PRESENT WOIK ; 

Present work is an attempt to reduce tlie complica- 
tions of three-dimensional analysis and yet simulate the 
steady state flow -field inside a single phase -shell and 
tube heat exchanger, A combination of volume porosity and 
distributed resistance is used in the simulation process. 

The effect of diffusion is also considered in the conserva- 
tion equations, A continuum approach is followed throughout 
the calculation domain. 
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The analysis is done iising two-dimensional conser- 
vation equations modified to represent the process in a 
shell and tube heat exchanger. The momentum and continuity 
equations are solved simultaneously using finite -element 
technique . The whole domain is divided into ei^t noded 
isoparametric elements , The Galerkin Weighted Residual 
method is used for obtaining the flow field and a three 
Gauss -point Numerical Integration is used for evaluating 
various integrals . Pressure f onnulation is used to solve 
the equations. The density of shell -side fluid was taken 
to be uniform throughout the domain, though the inclusion 
of their variations in the calculation procediare will present 
no difficulty , 

The objectives of the present work are : 

1 , To simulate the flow field inside single phase flow 
shell and tube heat exchangers . 

'2, to explore the possibility of using the above 

mentioned approach for effective modeling of heat 
exchangers . 
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CHAPTER-III 

FINITE ETJiMENT FORiyrULATIQN OF FLOW EQUATIONS 

3.1 PROBLEM DEFINITION : 

Figiire 3,1 shows a shell and tube heat exchanger 
with both the shell and tube sides arranged for a s ingle - 
pass flow. 

Shell -side fluid enters the heat exchanger through 
inlet (N-1 ) and comes out of the heat exchanger througli 
outlet (N-2) . The orientation of heat exchanger with respect 
to x,y and z axes is also shown in the same Figure. From 
the geometry of the shell and tube heat exchanger shown, 
one can see that the major directions along which the 
flxiid flows occur are X and Y directions . The flow in 
Z- direction is negligible. 

Further, the tube side fluid enteih the heat 
exchanger at a uniform temperature in the X-direction. 

The shell-side fliiid enters perpendicularly, also at a 
uniform temperature in the Y-direction. Hence the first 
row of tubes, which is at uniform temperature, will come 
in contact with shell -side fluid which is also at uniform 
temperature. Thus after passing the first row, the shell- 
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side fluid will be more or less at a uniform temperature . 
Hence it can be concluded without comitting a serious 
error idiat the shell -side fluid has temperature gradients 
only in the X-and Y-4irections , For the same reasons, the 
tube “Side fluid will also have temperature gradients in 
X and Y directions only. Hence the possibility of flow 
due to temperature gradient or natural convection currents 
is reduced in the Z -direction. 

For the reasons mentioned above, one can reduce 
the problem to two dimensions though with introduction of 
a slight error. For the tiibe side fluid the flow can be 
safely considered as unidirectional since the diameter of 
tube is generally small to allow the approximation. 


3.2 GOVERNING DIFFERENTIAL EQUATIONS : 


Figure 3,2 shows the diameteral section of a shell 
and tube heat exchanger in the X-Y plane. For simplicity, 
the inlet and outlet nozzle lengths are neglected, 
momentum and the continuity equations which govern the 
steady state flow of the shell -side fluid are given below. 
X-direction momentum equation ~ 


11 


^.3u 






« • • • 


(3.1) 
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Y-direction momentum equation - 




+ V ^ ^ ( d\ ^ s 


Ps^o 


(3.2) 


continuity equation - 


au . 3v _ ^ 
3x + 3y - ° 


(3.3) 


Where each term in the above equations is non- 
dimensionalized using the following parameters 


X = ,y 




= h ^ 


u = 


u. 


V = 


v_ 


and p = y 


(3.4) 


p„u 

o 


The left hand side of the momentum equations 
represents the convection terms and the second order terms 
on the right hand side represent the diffusion terms , The 
body forces appearing on the right hand side include (1 ) 
the weight of the fluid and (2) the distributed resistance 
in the x or y direction due to the presence of the tubes. 
It may be noted that the density appearing in the above 
equations is a fraction of the element volume filled by 
the shell -side fluid. This takes into account volume 
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porosity of the element. The porosity will, thus, have 
a value in the range of 0 to 1 . Throughout the analysis, 
the porosity is assmed to he a continuous parameter just 
as any other variable. In essence, the, concept of 
porosity essentially smears the effect of flow 
obstruction over the complete element volume. The location 
of the obstruction within the control vol\mie is not 
specified. The porosity, as used in the current analysis, 
is then a scalar quantity and does not vary with spatial 
coordinates , 

3.3 DISTRIBUTED RESISTANCE ; 


The following discussion presents a general 
description of the models used for the calculation of 
distributed resistances (Patankar, et .al., 1980) , 

As mentioned earlier, the effect of obstruction due 
to the presence of tubes is modeled using the concept of 
distributed resistance which represents the force of 
obstruction per unit volume of the shell side fluid. These 
resistances are present in the X and Y-directions only, 

3.3.1 X-Direction Resistance (R^) : 

This resistance accounts for the pressure drop due 
tc flow over the tube in ihie axial direction. This is 
calculated from (Patankar, 1980) . 
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( 5 - 5 ) 

3,2«2 Y-Direction Resistance i^y) ' 

This resistance accounts for the pressure drop due 
to flow perpendicular to the length of tube i,e. Hiis is 
a cross -flow resistance and is calculated uaing 
(Patankar, 1980) , 

fV ' -4 2^ .... (3.6) 

^ P 

PoV.t 

where G = (3,7) 

The resistances defined above are included in the 
body force terms in the conservation equations, along with 
tSie weight of the fluid which may act in any direction 
depending upon the orientation of the heat exchanger. 

3.4 FORMULATION FOR THE SHELL SIDE : 

The governing differential equations (3,1), (3.2) 
and (3,3) are solved simultaneously to obtain velocities 
inside the heat exchanger. 
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Figure 3,3 shows the calculation domain along with 
the boundary conditions for the shell side , No slip 
boundaries aire considered so that along the boundary, 
velocity in any direction is zero. Further, the fluid is 
assumed to enter the heat exchanger at a uniform velocity . 
The calculation domain is divided into a number of elements 
as shown. Hence u and v-velocities are known at any point 
on the boundary , Further the pressure at the inlet is 
known. 

The variables u, v and p are represented for a 
particular element using a polynomial which is function of 
the local coordinates for that element. These are nothing 
but the shape functions which are evaluated at Gauss points . 
Following an accepted practice of depicting the variation 
in pressure by shape functions of one order lower than those 
for defining the velocity distribution. 

n 

u = S N. u. 
i=1 ^ ^ 

n 

V = S N. V. 

i=i i i 

m 

p = E M. p. 
i=1 ^ ^ 


... (3.8a) 

. ... (3.8b) 


.,,.(3. 8c) 
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where n = 8 and m =4 

The requirement is fulfilled by utilising a 
rectangular isoparametric element with eight nodes where 
all the eight nodes are associated with velocity and only 
corner nodes are associated with pressure , Spatial 
co-ordinates within the element are defined in terms of 
co-ordinates of the eight nodes. The element is shown in 
Figiire 3,4. 

3,4,1 Shape Fmctions ; 

As stated earlier, shape functions are weighting 
functions, used to depict the variation of a variable in 
a particular domain. Once the calculation domain is divided 
into a number of elements, each element has certain number 
of nodes. The values of a particular variable on these 
nodes are used to define the variation of that variable 
inside the element (see equation (3.8)), 

Referring to Figure 3.4, the element shown has 
eight nodes, four corner nodes and four mid -side nodes. 

It may be noted that each side of the element has three 
nodes . Hence this element can accomodate a parabolic 
variation in the variable along any side of and across the 
element and is called a parabolic element. 





Y 


7 


6 



Sample 

Gausspoint 
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Local Co-ordinate System 
X-Y Global Co-ordinate System 


Figure 3-4 Eight Noded Rectangular 
Element with Nine Gauss 


Points 
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■ Shape function is a polynomial defined using a local 

co-ordinate system for the element. This is a normalised 
co-ordinate system, defined for each element in the 
domain. The shape functions used here are 

i) Corner nodes 

^ (1+ + n^n-1) .... (3.9a) 


ii) Mid -side nodes 

= 2 "l = 0 


;3.9b) 


Ni = I (1+ )<.^-■nh, 



.... (3,9c) 


Where and are local co-ordinates of i"*^ 

node , 

To relate local co-ordinate system with global 
co-ordinate system, again the same shape functions are 
used i.e, the variation of global co-ordinates inside an 
element is depicted b^? the same shape functions . Hence 


X = I N. 

.... (3.10a) 

y = I 

.... (3,10b) 
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^■fliere x and y refer to global co-ordinates of a 
point inside an element and x^^ and are global co-ordinates 
of i^^ node , It may be noted that the value of a shape 
function evaluated at a nodal point will be unity, since 
at a nodal point the value of the variable is already 
specified, Ihe element used here is isoparametric because 
same shape functions are used to define the variation in 
global co-03xiinates and the variable variation within an 
element. Once the shape functions are evaluated at a 
point inside an element, the values of any variable as well 
as global co-ordinates at that point can be determined. 

In present case (see Figure 3.4), there are nine Gauss 
points inside an element and each Gauss point will have 
eight shape functions . Hence there will be seventy two 
shape function values for each element. 

3.4,2 Variable Derivative : 

If a derivative of particular variable inside an 
element is required, it can readily be calculated using 
shape functions , Thus 


3u 


8 

= E 
i=1 


3N. 


• • • * 


(3.11) 
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Since nodal values of a variable are constant. 

The derivatives of shape fimctions with respect to 
global co-ordinates (see equation (3.11)) can be calculated 
using chain rule of partial differentiation. 




3x 

36 


3Nj_ 

W 



(3.12a) 


3N. ^ 

WT 3x ab 3y an 


.... (3.12b) 


which may be written as 


3N^/3e 


Sx/ae 

3y/3e 

3Nj^/3n 


3x/3n 

ay /an 


aN^/ae 


3N^/ax 


= J 


3N^/an ^ 


3N^/3n 


3N^/3x 

3Nj_/3y 


.... (3.13) 


Inverting J to find the global variation in shape 
functions, 


aN^/3x 

-1 

aNj_/3 6 

3N^/3y 

= J 

r 

3N^/an 


(3.14) 
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Where J is Jacobian and can be evaluated explicitly 
since the local variation in x and y can be defined. 


J= 


8 BN. 


8 3N, 


w ^ ifi ^ 


8 BN. 


^L^(gT,^)Xi 


8 BN. 

^.1 ^ Brr^ ^i 


.. ( 3 . 15 ) 


The inverse of Jacobian can be foimd out using 
standard matrix inversion technique. 



3 e 

dTl 

9x 


1 

!_ 

It 


3 6 

3y 

an 

1 

det J 

OX 

an 

3x 


.... ( 3 . 16 ) 


det J = i. IS - IS 


SX 

36 


( 3 . 17 ) 


It may be noted that 
dx dy = (det J) de d^i 


.... ( 3 . 18 ) 
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3,4,3 Algebraic Equations. 


Employing the Galerkin’s weighted residual approach, 
in which same approximating functions are used for the 
weighting and trial fimctions and using equation (3,8), 
equation (3,1) becomes 


ne 

E 

1 



n 

E 


Vk 


n 3N. n 

I 3? “ 3 + ^ Vk 



u 


3 


m 
+ E 


31^ 

3x 


51 - 


p u^ 
•^s o 


. n 3^N. 

^ ( S 

1 3x'^ 






n 

E 

1 


3^N. 


ay‘ 


u .) jr dA®= 0 


(3.19) 


Where the outer summation refers to each element in the 
domain and the inner siimmations to the appropriate nimber 
of nodes in an element. In the present case, n=8 and m=4. 

Invoking Green’s theorem, the order of the second 

order terms can be reduced leading to the expression, 

n 3^N. n 3^N. 

N. /“ E — u .+ E — ! 5 -a u. 7 dA® 

^ 1 3x^ ^ 1 3y^ 



f 


C® 


n 3N- 
N. E J 
^ 1 


371 


u. ds 


-] 


3N. n 3N. 




3Nj 

w 


n 

E 

1 


3y 


u._7dA' 


• • • • 


(3.20) 
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0 

Where C denotes the element siarface. VJhen these are 
summed over the adjacent element elements, the final 
contribution of the mass flux becomes zero lonless a boundary 
acts as a limit to the domain. 

Substituting equation (3,20) into equation (3,19) 


ne 

S 

1 


f 3N. 3N. 3M^ 

f j Wk ^ -3S Pi 


f L . 3N. 3N . 3N, 3N . ^ 7 


- ; 


Ps^o 


1 g 


( 3 . 21 ) 


where 


1 = 1 .... 4 
0=1 ..,,8 
K“1 ....S 

i = 1 .... No . of points . 


Ihe corresponding momentum equation in the Y- 
direction is obtained by simply interchanging x,y and u,v 
in equation (3,21). Substituting equation (3,8) into 
equation (3,3) yields, the continuity equation as 
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ne 

E 

1 



u. 


3N. 

* ^''3’ 


dA® =0 


.... (3.22) 


Where the 'weighting function is now taken as 
associated with four nodes only as far as the pressure is 
concerned. 

Ihe above equations are written for each element 
of the calculation domain and are solved sim'ultaneously to 
obtain the velocities in the domain. 

The assembled matrix equations take the form 


AX = F+B 


.... (3.23) 


where for the node i, the chosen form for X is 


Each coefficient in the matrix A has the form 
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^11 

^12 

S 3 

ne 

f 




^13 = f 

A® 

*^21 

*^22 

^23 



^31 

S 2 

S3 


1 




.... (3.24) 


where 


oN. BN. 

Sl= VA 33?^ 3 / 


3N. 3N. 3N. 3N. 


+ ‘ ( ^ __j. + __2: a > 

Re- '■ 3x 3x 3y 3y ' 


= 12 = "i 33r - 


°13 ■ '^21 “ "i 3x 


^99“ 


^23 =^i”Ty' ^31 = ° 


3M. 

^ 32 "% ^ 


C33= 


.... (3.25) 
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In the expression for and denote the 

velocities obtained from the previous iteration or the 
guessed velocities, since the equations to be solved are 
non-linear. In the present case, since all the boundary 
conditions are of essential type, the surface integral in 
equation (3,24) is not required and can be equated to zero. 

The first matrix on right hand side corresponds to 
the body forces and is written as 

dA® (3.26) 

f , L 

f =0, f, = N. .... (3.27) 

Ps^o 

■ The second matrix on the right hand side consists 
of the velocity gradient terms on the- boundary . In the 
present case, complete boundary is known and therefore 

this matrix is not required. It can be seen by examining 

- 9N - . ^ 

terms such as N. N, u, 3 that the interchange of i and 

1 K K 

3 results in a non— symmetric matrix. 


ne 

fl= s 


1 p 


where f^ = N- ^ , 


■1 
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3,4,4 Elemental Matrix : 

Referring to Figure 3,4, the element showi has 
ei^t nodes and it is already stated that corner nodes are 
associated with pressures and velocities where as midside 
nodes are associated with velocities only. Hence each 
corner node has three degrees of freedom (u-velocity, 
v-velocity and pressure) and a mid -side node has only two 
degrees of freedom (u-velocity and v-velocity) , dhus 
each element will have 20 degrees of freedom. The momentum 
and the continuity equations when wit ten for an element 
will give rise to a 20x20 matrix on left hand side as show 
below. 


1st Row 
2nd Row 

3rd Row 



I 


^20,1 ^20,2 ^20,3 ^20,4 


^1 , 5 ' * * • ^ 1,20 

^ 2,5 •*** ^ 2,20 

3,5 


® 20,5 ^ 20,20 

d 

(20x20) 
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V/here the coefficients are already defined in equation 
(3.24). 

In' above matrix 1st, 2nd and 3rd rows are obtained 
by writing X-direction momentum, continuity and Y-direction 
momentum equations respectively . The variable matrix for 
each element will be a column matrix of (20x1) as shown below. 


Pi 

^1 

u^ 

^2 

P3 

^3 


!_ 



(20x1 ) 


3,4,5 Global Matrix : 

The elemental matrix defined in section 3,4.4 is 
obtained for each element and all such matrices are 
assembled to obtain the global matrix. It may be noted that 
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while assembling the elemental matrices, some nodes will 
be common for more than one element. In such a case, the 
contributions from all the elements for that particular node 
are added and global matrix element for that node is 
calculated, Ihe global matrix will be a square NxN 
matrix of order N, where N is the total number of degrees 
of freedom for all the nodes in the calculation domain, 

3.4.6 Boundary Conditions : 

In the present case, velocities are known along the 
entire boundary and pressure is known at the inlet. To 
illustrate the incorporation of boundary conditions in 
global matrix, a (4x4) matrix is considered below. 

^11 ^12 ^13 ^14 ^ 

^21 ^22 ^23 ^24 ^2 = ^2 .... (3,26) 

^31 ^32 ^33 ^34 ^3 ^3 

^41 ^42 ^43 ^44__ [ "^ 4 ] L^4 _^ 

Assume that u^ is known and is equal to u^, then the 
incorporation of this boimdary condition in global matrix 
yields . 
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I 


— 



— 


~ - 


- 

1 


0 

^13 

^14 




F^ - a^ 2 

u 

0 

1 

0 

0 


^2 

= 

1 

u 


^31 

0 

^33 

^34 


^3 


Fj- 

t 

u 

^41 

0 

^43 

^44 


;^4 


V ^4 

1 


which can be reduced to 


^11 

a^3 

^14 

^31 

^33 

^34 

^41 

^43 

^44j 


— ^ 





F^- a^2 ^ 

^3 

= 

F3- a^3 u‘ 

U4 

— 

. 

F 4 - ^14 


(3.29) 


For a zero value at the boundary, equation (3,29) 
can be written as 


^11 ^13 ^14 





^31 ^33 ^34 


^ 3 : 

■ = 

^5 

^41 ^43 ^ 44 _ 



■ 

V '• 


3,4.7 Numerical Integration 

For evaluating the coefficients a^^^ in equation 

(3.24), the numerical integration is carried out using a 

3x3 Gauss point scheme. The sample gauss points are shQVhi 

in Fig\ire 3.4, The method used to carry out the numerical 

integration is as follows . 
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For each element the integrand, can he witten as 
- 1-1 +1 

I = y I F( e, 77) d£ d'n (3.31) 

-1 -1 - 

where e and 77 are local normalised co-ordinates. The 
integrand can he v;ritten as 

+1 +1 

I= y ( J F( e,T7) de )d7^ ....(3.32) 

-1 -1 - 


First carrying out the integration within the 
parentheses 


+1 

= / /■ S a. P( ^^,n)Jar) 


(3.33) 


The integrand is now a function of 77 only and 


hence 


m. 


"2 m 
I = L a. S' a- 
.1=1 ^ i=1 ^ 


1 - p(^i."i) 


.... (3.34) 


or 


j=m.^ m^ 

I = 2 E a. a. F( 77 ) 

3=1 1=1 J ^ ^ 


.... (3.34) 
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Here a 3x3 Gaussian integration is adopted hence 

= 3, where a^ and a. are the weighting factors for 
• th ^ , .th „ 

1 ana j Gauss points respectively. 

3.5 ITERATIVE PROCEDURE ; 

The set of non-linear simultaneous equations given 
by equation (3,23) is solved by an iterative process in 
which a simple convergence sequence and a method of 
"'^^^i^ble updating is employed. The procedure can be 
summarised as : 

i) Assume initial values of the variable u, p and v 

ii) Solve the matrix equations (3.23) to obtain 
new values of u, p and v at node points, 

iii) Evaluate the differences (u— u ) /u, (v— v*) /v and 
(p-p)/pl If thes-e are within the specified 
tolerance at all points then the calculation is 
complete, 

iv) If the differences calculated above are not 

within the required tolerance then using suitable 
under-relaxation determine the new values of 
u, p and V. 

The steps (i) to (iv) are repeated until convergence 


is obtained. 
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CHAPTER- IV 


PROGRAMS AND PROCEDURES 


4.1 SUBROUTINES : 

The program is made up of following subroutines. 

1 . DIMENS . 

A number of vectors and arrays are utilized during 
house-keeping process involved in the formiiLation of the 
matrix to be solved. The dimension required for each of 
these vectors or arrays is set in subroutine DIMENS, This 
permits a form of dynamic dimensioning to be \ised. The 
program can, therefore, be increased or decreased in size 
to suit a particular grid size,. 

2 . DINPUT. 

The data required are read in subroutine DINPUT, 

This includes physical properties of the shell -side fluid 
and boundary conditions. As the data are given to two 
other subroutines D1AGN1 and DIAGN2, they are csilled in to 
check that the data which are given comply with the 
required checks . 

3 . DIAGN1 . 

This subroutine checks the overall control parameter 

which governs the nimiber of nodes, number of elements and 

the boundary conditions for the particular problem under 
investigation , 
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4. DIAGN2. 

Various checks are incorporated to ensure that the 
geometric data obey simple criteria. If any check fails, 
the data which are not read to that point will be read and 
printed , 

5. DRIVES. 

A subroutine in which the constants necessary for 
numerical integration of the relevant equations are defined. 
Subroutines SHAPES, SHAPE4 and DJACOB are called from DRIVES . 

6. SHAPES. 

Ihe subroutine SHAPES calculates the shape functions, , 
required for defining the variation of the independent 
variable over an element which contains eight nodes on the 
element boundary, 

7. SHAPE4. 

This subroutine calculates the shape functions for 
a four noded element. Two shape function subroutines are 
used in the program since the pressure and velocity in the 
Navier -Stokes equations here are associated with four and 
eight nodes of an element respectively. 

8. DJACCB. 

In this subroutine, the first order differentials 
of the shape functions, with respect to the chosen co-ordinate 
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system, are calculated. These are required v^en terms such 
„ 3u 3u , 

9x ’ '5y calculated. For instance, having 

defined . 


Z”^'] » ^2 , , Nq_J7 


U2 

^3 


then first order differentials can he written in the form 


3U ^ 3Ng 

3x 3x * 3x * 3x ’*“** dx 




# 


u 


8 


since N is a function of spatial co-ordinates x,y. 

9. ITERATi 

'This is the main subroutine which calls the 
necessary subroutines for the evaluation and solution of 
the necessary matrix equations, transformation of boundary 



37 


conditions into a suitable form for incorporation into the 
matrix equations, output of results and check on the 
convergence tolerances when an iterative procedure is 
required, 

10. FRONTS. 

This subroutine formulates the global matrix, imposes 
the .boundary conditions and solves the resulting matrix 
system of equations using a non-symmetric Frontal meldiod. 

The required global matrix is assembled from contributions 
of element matrices as set up in the subroutine MATRIX, 

11. MATRIX, 

This subroutine forms the necessary relevant fluid 
matrix and corresponding right hand side vector, 

12. TOLREL. 

In this subroutine a check is made to ensure that the 
tolerance limits are met. If these are not met, updated 
values of the required coefficients are evaluated and the 
solution procedure is repeated, 

13. MESGENi 

This subroutine is used for the calculation of space 
co-ordinates for each node in the mesh. The input for this 
is given in the form of space intervals in x and y directions. 



Master Fluid 



General Program Structure 



















Fig.4.2 Flow-Chart 
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4.2 GENERAL STRUCTURE OF THE PROGRAM : 

The general structure of the program alongwith 
various subroutines is shown in Figure 4,1 . The flow-chart 
is shown in Figure 4.2. The working of the program 
consists of the following steps, 

2* The subroutine DIMENS is called from the main 
program MASTER FLUID. This will set the dimensions of 
various arrays such as nxmiber of points, number of degrees 
of freedom in the program, 

2. Subroutine DINPUT is called from the main program. 
This subroutine will collect information such as initial and 
boundary conditions of the problem. In DINPUT the 
program goes through the following steps . 

a. Subroutine MESGEN will be called from DINPUT 
for calculation of nodal co-ordinates and mesh- 
generation, 

b. Subroutines DIAGN1 and DIAGN2 are then called to 
check the data given throiogh input, satisfy certain 
rules such as number of boundary conditions etc. 

3 , Subroutine DRIVES is called from the MATRIX 

This subroutine will calculate shape functions. Jacobians 
and global derivatives of the shape functions. The various 

steps in this are: 
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a. Subroutine SHAPES and SHAEE4 will be called to 
calculate the shape functions for each element, 

b » Subroutine DTACOB will be called to calculate the 
derivatives of the shape functions with respect to, 
the local co-ordinates, 

c. Global derivaties of tiie above shape functions will 
be calculated, 

d. Subroutine SHAPE4 will be called to calculate 

shape functions for 4 noded element, 

e. Global derivatives of above shape functions will 
be calculated, 

4. Subroutine ITERAT will be called to obtain the 

solution, The procedure is given below. 

a. Subroutine FRONIS will be called to solve the 

governing equations. Here MATRIX subroutine is 
called for each element and these elemental matrices 
are assembled in subroutine FRONTS. The assembled 

matrix equations are solved and new values ofthe 
v3jriSLbl6S siTG • 

b. After calculating new values of variables, Subroutnn' 

TOLiI^L will be called to compare the difference 

and the old values of variables. If the 


between new 
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solution is within the required tolerance, control 
goes to subroutine WRITER to write the converged 
solution. 

If the solution obtained is not within the requiri^^ 
tolerance, new values of variables are relaxed using some 
relaxation factor and are taken as the new guesses for 
fresh iteration. The fresh iteration starts with a call 
to subroutine FRONTS . 

4,3 SOLUTION OF ALGEBRAIC EQUATIONS : 

As discussed in Chapter 3 the governing differential 
equations are converted into algebraic equations. 

First of all, the calciolation domain was divided 
into a number of elements and in each eleL-ent certain 
number of points were selected as representatives of the 
pressiare and velocity field inside the element. Here, 
eight points are chosen to represent velocity field and 
four points are chosen to represent pressure field inside th 

element . 

When the governing differential equations are 
descretized for each element, a matrix of order 20 will 

solution will give values of variables 


be obtained whose 
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at the nodes , When such elemental matrices are assembled, 
the resulting global matrix in the present case becomes 
unsymmetric. Straight forward assembly of such elemental 
matrices requires a huge core memory in the computer 
when the number of elements increase* Hence, in 'the present 
case Frontal method is used to solve the matrix equations. 
The choice of the method has a significant bearing on computer 
storage requirement and execution time. 

The overall solution technique consists of the 
following steps, 

(i) Formation of elemental matrices (in Subroutine 
MATRIX) 

(ii) Assembling into a global matrix. 

(iii) Introduction of boundary conditions. 

(iv) Reduction of the global fluid matrix \ising a 
Gaussian elimination technique, 

(v) A back -substitution process. 

First step in above sequence is already discxissed. 

When the assembly of elemental matrices starts, then before 
the global matrix is completely assembled, certain 
equations can be eliminated using Gaussian elimination 
technique. This principle is used in the Frontal method. 
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As soon as all contributions from all -Uae relevant 
elements to a particular nodal point are assembled then 
corresponding variables associated with that node can be 
eliminated. The complete matrix, is therefore, never 
assembled and the eliminated equations can be stored on a 
disc. The equations held in the core, with corresponding 
nodes and variables, are termed as the front and number 
of unknovn variables in front are termed as front width. 
Thus the order of the global matrix solved at the end will 
not be more than the specified front width. After the 
global matrix is solved, the eliminated equations stored 
on the disc can be solved using this solution. 
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cRkpim- V 

RESULTS, AND DISCUSSION 

5.1 GENERAL GUIDELINES : 

As an example, the calculation domain divided into 
64 elements is shown in Figure 3,3. The mesh is uniform 
in the Y-direction and non-uniform in the X-direction, 
Further, while dividing the mesh, care was taken to ensure 
that the aspect ratio of element does not exceed 8. This 
affects the accuracy of the solution. 

Further, the ratio of the number of unknown 
velocities to the number of unknown pressures was atleast 2, 
If the ratio became less than two, it resiT.ted in an ill- 
conditioned matrix. This was because the original equations 
contained second order derivatives of velocity and first 
order derivatives of pressure . To achieve this x^tio, the 
element chosen was such that it had four pressure points 
(only comer points) and eight velocity points (comer 
points as well as mid-side points). Thus the mesh was 
designed to maintain this ratio (i.e, atleast 2 as stated 
above) . 
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In the present case, the pressure could be specified 
only at the inlet and hence the mesh was finer here, so that 
larger number of pressures could be specified. The inlet 
velocity profile was assumed parabolic so that velocity was 
maximum at the centre of inlet and zero at the ends , 

The results for various values of the parameters of 
a shell and tube heat exchanger are shown in Figures 5,1 
through 5,7, Some explanation is given in each case. However, 
as no data were available to compare the results, the 
explanations are only qualitative . The results are used for 
validating the model rather than for predicting the performance. 
For all the cases presented here following data were used. 

Tube diameter ip mm 

Tube pitch (square arrangement) 25 mm 

Ratio of Shell length to shell 2 

diameter . 

Ratio of inlet pipe diameter 0,1 

to shell diameter. 

Heat exchanger inlets and outlets were taken to be 
of the same size, though any combination could be tried. 

The length of the heat exchanger and the center line 
velocity in the inlet pipe were used to non dimens ionalize 
the equations. Hence the velocities calculated were non- 
dimensional , 
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For tho heat exchanger tuho-bank, "the friction 
factors are given below(Salisbury, 1895). 

(i) For flov/ along the tube bank, f= 0.02 

(ii) For flow across the tube -bank, fr= 0,24. 

The friction factors used here are applicable only 
when the tubes are arranged in a square pitch. It is to be 
noted that for calculating ihe pressure drops the equivalent 
diameter used is tho hydraulic diameter which is defined as 
four times tho free liquid -passage area divided by the liquid 
touched perimeter of tho tubes . 

Tho effect of changing tho orientation of tho heat 
exchanger and the Reynold's number is discussed in tho following 
sections . 

5.2 CASE~I : 

Hero, ' the fluid enters and loaves tho heat exchanger 
in the direction of gravity i,e. gravity is in tho positive 
Y-diroction. Tho results are shown in Figures 5,1 through 
5.5 . Tho value of Reynold’s number was taken to be 1000, 

Figure 5.1 shows velocity distributions along tho 
lino A-A (see Figure 3.3), The u-velocity is maximum after 
some distance from tho inlet and then keeps on decreasing 
\jji-til ii; becomes zero at the wall. The v— velocity is maximum 



at thG inlot and then suddenly drops as the liquid- expands 
after entering the heat exchanger. Then after a sILigii't 
increase it also keeps on decreasing until it becomes zero at 
the wall . 

Figure 5^2 shows velocity distributions alc^ng the line 
B-B (see Figure 3,3), This being the transverse ;^ection 
taken at the center of heat exchanger, u-velocity is maximuia 
at the center and it reduces to zero at the walls ^ Further 
because of the position of inlet, the u-velocity ext a parti- 
cular point on B-B downstream of the point Ogis g 3 r*eater 
than that at the corresponding point upstream of "tiiG point Og. 
The v-velocity becomes maximum after some distance from the 
wall, then it reduces and remains almost constant in the 
central portion of the section. Near the other eXii of the 
wall ’ { RS ) it attains a slightly negative value loef ore 
reaching zero at the wall. 

Figure 5,3 shows velocity distributions a3-ong the 
line C-C which is a longitudinal section near the inlet of the 
heat exchanger (see Figure 3,3). After entering "th^ heat 
exchanger, the liquid flows in all directions and- hence 
u-velocity becomes negative for some portion alors-g C-C, It 
then increases and remains almost constant for s<a>Eie length of 
the heat exchanger. It reaches zero at the wall after a slight 
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incroasG, This fluctuation is because of eddies near the 
wall. For v-velocity two distinct peaks are seen at equal 
distances on both sides of the center line of inlet. This • 
is because of the sudden esqjansion of idie liquid which makes 
it travel side -ways leaving the centre line. The right hand 
side peak is slightly larger 'vdii.ch may be due to the fact 
that the outlet is on the right hand side of the section hence 
tendency of the liqx:iid to travel more to ihe right. After 
second peak, v-velocity decreases and changes its direction 
several times before reaching a zero value due to the "eddies 
prevailing near -the wall , ' . , 

Figure 5,4 shows velocity distribution along the ' 
line D-i) (see Figure 3,5) which is a longitudinal section ' 
passing through the center of the heat exchanger, [Hie 
velocity distribution is symmetric about the center, Ihe 
u-velocity attains a maximijm value near the center while the 
v-velocity reaches a minjmnTn there. 

Figure 5,5 is a plot of stream lines which shows an 
expected pattern. Stream linos having very low or high values 
of stream feinctions are closer to the walls of the heat 
exchanger. The eddies (not shown) are predominant near the 
wall PQ(i.o, wall opposite to the outlet of the heat exchanger 
Further, there are almost no eddies near the outlet comer S , 
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The corner R opposite to the inlet wall has small eddies 
present (not shown) , 

5.3 CASE -II . 

For this case, the orientation of heat exchanger was 
changed with respect to the direction of gravity i,e. the 
gravity acts in the positive X-<iirection , Figure 5,6 shows 
the nature of stream lines in this case for the same value of 
Reynold's number of 1000. It may be noted that there is 
a sli^t shift of stream lines in the X-direction as 
compared to the previous case. But it may be noted that "there 
is a very little effect on the flow field because of change 
in the. -direction of gravity. Few runs were taken by changing 
the direction of gravi"fcy in negative Y and X-directions . 

It was noted that if the direction of gra-^rity was changed 
from positive Y-direction to negative Y-direction» "the 
flow-field showed almost insignificant change. On the 
other hand, >hen "the direction of gravi"ty was changed from 
positive X-direction to negative X-direction, there was 
a noticeable change in the flow-field. This may be explained 
as follows. The Y-force consists of resistance due to cross- 
flow over the tubes and gravi"ty force if present in that 
direction. Now cross -flow resistance being much larger in 
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value dominates the gravity force. Hence even if the 
direction of gravity is changed, it has no significant 
effect on the flow -field, Mhereas the resistance for axial 
flow over the tubes is much smaller in magnitude compared 
to the cross -flow resistance and it is comparable with the 
force of gravity, v^ich explains the change in the flow -field 
with change in the direction of gravity frcan positive to 
negative X-direction, 

5.4 CASE -III. 

Here it was assumed that the gravity was in the 
positive Y-direction, but the value of Reynold’s number was 
changed to 2000, 

Figure 5,7 shows the velocity vectors in the 
calculation domain of the heat exchanger. It was noted that 
the nature of flow-field remained almost same except slight 
increase in the eddies inside , In this case the effect of 
gravity on the flow field was even more insignificant 
compared to the previovis cases. It can be seen that as the 
Reynold’s number increases the effect of gravity becomes 
even move insignificant. 


'P 
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5.5 CONCLUSIONS : 


A FEN program was developed to solve the governing 
equations given in Chapter 3. The solution yielded the 
velocity distribution in a shell and tube heat exchanger 
without phase change. The computed res\alts showed that 
for a given aspect ratio, the solution diverged for a parti- 
cular value of Reynold* s number. As the aspect ratio increased 
the limit of Reynold's number beyond vhich divergence 
occtired was also increased. The program can be modified 
for use at high Reynold's number by incorporating eddy 
viscosity using the method given by Sha et, al. ^'1982^7. 

Despite the difficulty of getting convergence at 
high Reynold's number as mentioned above, the results obtained 
in this thesis can be considered to be pretty realistic as 
the aspect ratio in general for shell and tube heat exchangers 
is much greater than two. The range is three to fifteen 
but a value of eight is generally used, Bhaskare, 1986_7 
which will make it possible to get the converged solution 
at higher Reynold's number. 

The objective of this work was to develop a code 
for flow-modeling of a shell and tube heat exchanger without 
phase change . Since this is the first attempt of two 
dimensional flow-modeling in a shell and tube heat exchanger 
using FEM,only the simpler cases of laminar flow at lov/ 
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Reynold’s number are dealt with in this thesis. There is 
a large scope for further work on this project as well as 
for refinements in the computation. The results are used 
only to validate the model as no data are available for 
ccMparison. 

5.6 RECOMMENDATIONS FOR FUTURE WORK : 

The code in the present version is limited to the 
two dimensional analysis of a shell and tube heat exchanger 
without phase change. The code developed has a flexibility 
to accomodate various changes easily. Since the method 
used is FEM, it can be extended easily to handle other 
geometries as well. For example, flow modeling for a_ plate 
type or spiral type heat exchanger can be easily done 
because of the fact that the element used in the code is 
parabolic, which can fit irregular boundaries. Further, 
the code can handle any type of boundary conditions i.e. 
variable specified boundary or the gradient specified 
boundary , 

With little effort, energy equation can be added in 
the present code to use it for analysing thermal performance 
of a heat exchanger. The elemental matrix will be 28x28 
in that case, with four degrees of freedom at the corner 
nodes and three degrees of freedom at the mid -side nodes. 
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For flow -modeling at very high values of Reynold's 
nijunher say above 5000, a suitable turbulent viscosity model 
can be added to the present code which will update the 
viscosities at eyery node depending on the velocity at that 
node , 

Further, the present code can be extended to 
three dimensional situations. 

In case I and II, for a tolerance of + 3 '/. , 12 iteratic 
were required and the CPU time was 8 minutes. In case III, for 
the same tolerance limit 15 iterations were required and the 
CPU time was 10 minutes. 
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